#ifndef	__MathLib_h__
#define	__MathLib_h__

#if _MSC_VER > 1000
#	pragma once
#endif	// _MSC_VER > 1000

#ifdef _WIN32
#	pragma warning(disable: 4244)
#	pragma warning(disable: 4267)
#	pragma warning(disable: 4305)
#endif

#include "inccommon.h"

namespace rtMath
{

const float	EPSILON		= 1e-6f;
const float	PI			= 3.14159265358979323846f;
const float	DBL_PI		= 6.28318530717958647692f;
const float	HALF_PI		= 1.57079632679489661923f;
const float	SQRT_PI		= 0.785398163397448309615f;
const float	INV_PI		= 1.0f / PI;
const float	DEG2RAD		= PI / 180.0f;
const float	RAD2DEG		= 180.0f / PI;
const float	MAX_FLOAT	= std::numeric_limits<float>::max();
const float	MIN_FLOAT	= std::numeric_limits<float>::min();

inline float rsqrt_f(float x)	// faster version
{
	union f2i {
		float	f;
		int32	i;
	} tmp;
	
	const float x2 = x * 0.5f;
	tmp.f = x;
	tmp.i = 0x5f3759df - (tmp.i >> 1);
	tmp.f = tmp.f * (1.5f - (x2 * tmp.f * tmp.f));
	return tmp.f;
}

inline float rsqrt(float x)
{
	union f2i {
		float	f;
		int32	i;
	} tmp;
	
	const float x2 = x * 0.5f;
	tmp.f = x;
	tmp.i = 0x5f3759df - (tmp.i >> 1);
	tmp.f = tmp.f * (1.5f - (x2 * tmp.f * tmp.f));
	tmp.f = tmp.f * (1.5f - (x2 * tmp.f * tmp.f));	// iOrange - I've added second Newton-Rapson iteration to the original code
	return tmp.f;									//			 for better result precision
}

template <typename T>
inline T clamp(const T& val, const T& _min, const T& _max)
{
	return (val < _max ? (val > _min ? val : _min) : _max);
}
template <typename T>
inline T sign(const T& val)
{
	return (val < T(0)) ? T(-1) : T(1);
}
inline int nextpow2(int n)
{
	n = n - 1;
	n = n | (n >> 1);
	n = n | (n >> 2);
	n = n | (n >> 4);
	n = n | (n >> 8);
	n = n | (n >> 16);
	n = n + 1;
	return n;
}


class vec2
{
public:
	vec2() : x(0), y(0) { }
	vec2( const vec2 &v ) : x(v.x), y(v.y) { }
	vec2( float _x, float _y ) : x(_x), y(_y) { }
	
	explicit vec2( const float* v ) : x(v[0]), y(v[1]) { }
	
	inline bool operator == ( const vec2& v ) const { return fabs(x - v.x) < EPSILON && fabs(y - v.y) < EPSILON; }
	inline bool operator != ( const vec2& v ) const { return !(*this == v); }
	
	inline const vec2 operator * ( float f ) const { return vec2(x * f,y * f); }
	inline const vec2 operator / ( float f ) const { return vec2(x / f,y / f); }
	inline const vec2 operator + ( const vec2 &v ) const { return vec2(x + v.x,y + v.y); }
	inline const vec2 operator - ( void ) const { return vec2(-x,-y); }
	inline const vec2 operator - ( const vec2& v ) const { return vec2(x - v.x,y - v.y); }
	
	inline vec2& operator *= ( float f ) { return *this = *this * f; }
	inline vec2& operator /= ( float f ) { return *this = *this / f; }
	inline vec2& operator += ( const vec2& v ) { return *this = *this + v; }
	inline vec2& operator -= ( const vec2& v ) { return *this = *this - v; }
	
	inline float operator * (const vec2 &v) const { return x * v.x + y * v.y; }
	
	inline operator float*() { return v; }
	inline operator const float*() const { return v; }
	
	inline float& operator[]( int i ) { return v[i]; }
	inline const float operator[]( int i ) const { return v[i]; }
	
	inline float length( void ) const { return sqrt(x * x + y * y); }
	inline vec2& normalize( void )
	{
		const float inv = rsqrt(x * x + y * y);
		x *= inv;
		y *= inv;
		return *this;
	}
	
	inline void saturate( void )
	{
		if(x < 0.0f) x = 0.0f;
		else if(x > 1.0f) x = 1.0f;
		if(y < 0.0f) y = 0.0f;
		else if(y > 1.0f) y = 1.0f;
	}
	
	inline void clamp( const vec2& min, const vec2& max )
	{
		if		( x < min.x ) x = min.x;
		else if	( x > max.x ) x = max.x;
		if		( y < min.y ) y = min.y;
		else if	( y > max.y ) y = max.y;
	}
	
	union
	{
		struct { float x, y; };
		struct { float s, t; };
		float v[2];
	};
};

class vec3
{
public:
	vec3()
#ifndef	RT_USE_SIMD_SSE
		: x(0), y(0), z(0)
#endif	// RT_USE_SIMD_SSE
	{
#ifdef	RT_USE_SIMD_SSE
		_sse = _mm_setzero_ps();
#endif	// RT_USE_SIMD_SSE
	}
	vec3(const vec3& v)
#ifndef	RT_USE_SIMD_SSE
		: x(v.x)
		, y(v.y)
		, z(v.z)
#endif	// RT_USE_SIMD_SSE
	{
#ifdef	RT_USE_SIMD_SSE
		_sse = v._sse;
#endif	// RT_USE_SIMD_SSE
	}
	vec3(float _x, float _y, float _z)
#ifndef	RT_USE_SIMD_SSE
		: x(_x)
		, y(_y)
		, z(_z)
#endif	// RT_USE_SIMD_SSE
	{
#ifdef	RT_USE_SIMD_SSE
		_sse = _mm_set_ps(_x, _y, _z, 0.0f);
#endif	// RT_USE_SIMD_SSE
	}

	explicit vec3(const vec2& v)
#ifndef	RT_USE_SIMD_SSE
		: x(v.x)
		, y(v.y)
		, z(0)
#endif	// RT_USE_SIMD_SSE
	{
#ifdef	RT_USE_SIMD_SSE
		_sse = _mm_set_ps(v.x, v.y, 0.0f, 0.0f);
#endif	// RT_USE_SIMD_SSE
	}
	explicit vec3(const float* v)
		: x(v[0])
		, y(v[1])
		, z(v[2])
	{ }

	inline bool operator == (const vec3& v) const { return fabs(x - v.x) < EPSILON && fabs(y - v.y) < EPSILON && fabs(z - v.z) < EPSILON; }
	inline bool operator != (const vec3& v) const { return !(*this == v); }
	inline bool operator < (const vec3& v) const { return x < v.x && y < v.y && z < v.z; }
	inline bool operator > (const vec3& v) const { return x > v.x && y > v.y && z > v.z; }
	inline bool operator <= (const vec3& v) const { return x <= v.x && y <= v.y && z <= v.z; }
	inline bool operator >= (const vec3& v) const { return x >= v.x && y >= v.y && z >= v.z; }

	inline const vec3 operator * ( float f ) const { return vec3(x * f,y * f,z * f); }
	inline const vec3 operator / ( float f ) const { return vec3(x / f,y / f,z / f); }
	inline const vec3 operator + ( const vec3& v ) const { return vec3(x + v.x, y + v.y, z + v.z); }
	inline const vec3 operator - () const { return vec3(-x,-y,-z); }
	inline const vec3 operator - ( const vec3& v ) const { return vec3(x - v.x, y - v.y, z - v.z); }
	inline const vec3 operator - ( float f ) const { return vec3(x - f, y - f, z - f); }

	inline vec3& operator *= ( float f ) { return *this = *this * f; }
	inline vec3& operator /= ( float f ) { return *this = *this / f; }
	inline vec3& operator += ( const vec3& v ) { return *this = *this + v; }
	inline vec3& operator -= ( const vec3& v ) { return *this = *this - v; }

	inline vec3 operator * (const vec3& v) const { return vec3(x * v.x, y * v.y, z * v.z); }

	inline operator float*() { return v; }
	inline operator const float*() const { return v; }

	inline float& operator[](int i) { return v[i]; }
	inline const float operator[](int i) const { return v[i]; }

	inline vec3& set(float _x, float _y, float _z) { x = _x; y = _y; _z = z; return *this; }
	inline float length( void ) const { return sqrt(x * x + y * y + z * z); }
	inline float lengthSq( void ) const { return x * x + y * y + z * z; }
	inline void normalize( void )
	{
		float inv = 1.0f / sqrtf(x * x + y * y + z * z);
		x *= inv;
		y *= inv;
		z *= inv;
	}

	inline void saturate( void )
	{
		if		( x < 0.0f ) x = 0.0f;
		else if	( x > 1.0f ) x = 1.0f;
		if		( y < 0.0f ) y = 0.0f;
		else if	( y > 1.0f ) y = 1.0f;
		if		( z < 0.0f ) z = 0.0f;
		else if	( z > 1.0f ) z = 1.0f;
	}

	inline void clamp( const vec3& min, const vec3& max )
	{
		if		( x < min.x ) x = min.x;
		else if	( x > max.x ) x = max.x;
		if		( y < min.y ) y = min.y;
		else if	( y > max.y ) y = max.y;
		if		( z < min.z ) z = min.z;
		else if	( z > max.z ) z = max.z;
	}

	union
	{
#ifndef RT_USE_SIMD_SSE
		struct { float x, y, z; };
		struct { float r, g, b; };
		float v[3];
#else // RT_USE_SIMD_SSE
		struct { float x,y,z,_padd; };
		float v[4];
		__m128 _sse;
#endif // RT_USE_SIMD_SSE
	};
};

inline float dot(const vec3& v1,const vec3& v2)
{
	return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
}
inline vec3 cross(const vec3& v1,const vec3& v2)
{
	return vec3(v1.y * v2.z - v1.z * v2.y,
				v1.z * v2.x - v1.x * v2.z,
				v1.x * v2.y - v1.y * v2.x);
}
inline vec3 normalize(const vec3& v)
{
	float inv = 1.0f / sqrtf(v.x * v.x + v.y * v.y + v.z * v.z);
	return vec3( v.x * inv, v.y * inv, v.z * inv );
}
inline vec3 reflect(const vec3& dir, const vec3& normal)
{
	vec3 temp = normal * dot(dir, normal) * 2.0f;
	return normalize(dir - temp);
}
inline vec3 refract(const vec3& I, const vec3& N, const float eta)
{
	return normalize((N + I) * eta - N);
}
inline float distance(const vec3& a, const vec3& b)
{
	vec3 c = a - b;
	return c.length();
}
inline float distance_square(const vec3& a, const vec3& b)
{
	vec3 c = a - b;
	return c.lengthSq();
}
inline vec3 transform_by_normal(const vec3& normal, const vec3& v)
{
	float sinV = -v.y;
	float cosV = sqrtf(1.0f - sinV * sinV);

	float cosU, sinU;

	if (cosV == 0.0f)
	{
		sinU = 0.0f;
		cosU = 1.0f;
	}
	else
	{
		float invCosV = 1.0f / cosV;
		sinU = v.x * invCosV;
		cosU = v.z * invCosV;
	}

	return vec3((normal.y * sinV + normal.z * cosV) * sinU + normal.x * cosU,
				 normal.y * cosV - normal.z * sinV,
				(normal.y * sinV + normal.z * cosV) * cosU - normal.x * sinU);
}


class vec4
{
public:
	vec4() : x(0), y(0), z(0), w(1) { }
	vec4(const vec4& v) : x(v.x), y(v.y), z(v.z), w(v.w) { }
	vec4(const vec3& v, float _w) : x(v.x), y(v.y), z(v.z), w(_w) { }
	vec4(float _x, float _y, float _z, float _w) : x(_x), y(_y), z(_z), w(_w) { }
	
	explicit vec4(const vec3& v) : x(v.x), y(v.y), z(v.z), w(1) { }
	explicit vec4(const float* v) : x(v[0]), y(v[1]), z(v[2]), w(v[3]) { }
	
	inline bool operator == (const vec4& v) const { return fabs(x - v.x) < EPSILON && fabs(y - v.y) < EPSILON && fabs(z - v.z) < EPSILON && fabs(w - v.w) < EPSILON; }
	inline bool operator != (const vec4& v) const { return !(*this == v); }
	
	inline const vec4 operator * (float f) const { return vec4(x * f,y * f,z * f,w * f); }
	inline const vec4 operator / (float f) const { return vec4(x / f,y / f,z / f,w / f); }
	inline const vec4 operator + (const vec4& v ) const { return vec4(x + v.x,y + v.y,z + v.z,w + v.w); }
	inline const vec4 operator - () const { return vec4(-x,-y,-z,-w); }
	inline const vec4 operator - (const vec4& v) const { return vec4(x - v.x,y - v.y,z - v.z,w - v.w); }
	
	inline vec4& operator *= (float f) { return *this = *this * f; }
	inline vec4& operator /= (float f) { return *this = *this / f; }
	inline vec4& operator += (const vec4& v) { return *this = *this + v; }
	inline vec4& operator -= (const vec4& v) { return *this = *this - v; }
	
	inline float operator * ( const vec4& v ) const { return x * v.x + y * v.y + z * v.z + w * v.w; }
	
	inline operator float*() { return v; }
	inline operator const float*() const { return v; }
	
	inline float& operator[](int i) { return v[i]; }
	inline const float operator[](int i) const { return v[i]; }
	
	inline vec4& set(float _x, float _y, float _z, float _w) { x = _x; y = _y; z = _z; w = _w; return *this; }
	inline float length(void) const { return sqrt(x * x + y * y + z * z + w * w); }
	inline void normalize(void)
	{
		float inv = 1.0f / sqrtf(x * x + y * y + z * z + w * w);
		x *= inv;
		y *= inv;
		z *= inv;
		w *= inv;
	}
	
	inline void saturate(void)
	{
		if		(x < 0.0f) x = 0.0f;
		else if	(x > 1.0f) x = 1.0f;
		if		(y < 0.0f) y = 0.0f;
		else if	(y > 1.0f) y = 1.0f;
		if		(z < 0.0f) z = 0.0f;
		else if	(z > 1.0f) z = 1.0f;
		if		(w < 0.0f) w = 0.0f;
		else if	(w > 1.0f) w = 1.0f;
	}
	
	inline void clamp(const vec4& min, const vec4& max)
	{
		if		(x < min.x) x = min.x;
		else if	(x > max.x) x = max.x;
		if		(y < min.y) y = min.y;
		else if	(y > max.y) y = max.y;
		if		(z < min.z) z = min.z;
		else if	(z > max.z) z = max.z;
		if		(w < min.w) w = min.w;
		else if	(w > max.w) w = max.w;
	}
	
	union
	{
		struct { float x, y, z, w; };
		struct { float r, g, b, a; };
		float v[4];
#ifdef RT_USE_SIMD_SSE
		__m128 _sse;
#endif // RT_USE_SIMD_SSE
	};
};

inline float dot(const vec4& v1, const vec4& v2)
{
	return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z + v1.w * v2.w;
}
inline uint32 vec4_to_uint32(const vec4& v)
{
	uint32 r = static_cast<int>(v.r * 255.0f);
	uint32 g = static_cast<int>(v.g * 255.0f);
	uint32 b = static_cast<int>(v.b * 255.0f);
	uint32 a = static_cast<int>(v.a * 255.0f);

	return (r << 24) | (g << 16) | (b << 8) | a;
}
inline vec4 uint32_to_vec4(uint32 c)
{
	static const float _1by255 = 1.0f / 255.0f;
	return vec4( (c & 0x000000ff)        * _1by255,
				((c & 0x0000ff00) >>  8) * _1by255,
				((c & 0x00ff0000) >> 16) * _1by255,
				((c & 0xff000000) >> 24) * _1by255);
}


class Ray
{
public:
	Ray()
		: lightStrength(1)
	{}
	Ray(const vec3& _origin, const vec3& _dir)
		: lightStrength(1)
		, origin(_origin)
	{
		SetDir(_dir);
	}
	Ray(const Ray& _other)
	{
		this->operator =(_other);
	}

	Ray& operator = (const Ray& other)
	{
		origin = other.origin;
		dir = other.dir;
		dir_inv = other.dir_inv;
		posneg[0] = other.posneg[0];
		posneg[1] = other.posneg[1];
		posneg[2] = other.posneg[2];
		return (*this);
	}

private:
	void inline SetDir(const vec3& _dir)
	{
		dir = normalize(_dir);
		dir_inv = vec3(1.0f / dir.x, 1.0f / dir.y, 1.0f / dir.z );
		posneg[0] = (dir.x > 0 ? 0 : 1);
		posneg[1] = (dir.y > 0 ? 0 : 1);
		posneg[2] = (dir.z > 0 ? 0 : 1);
	}

public:
	vec3	origin;
	vec3	dir;
	vec3	dir_inv;
	int		posneg[3];
	float	lightStrength;
};

template<typename T>
inline T SolveBarycentric(const T& t0, const T& t1, const T& t2, float u, float v)
{
	return t0 * (1 - u - v) + t1 * u + t2 * v;
}

// iOrange - angle in degrees
static void RotateVector(vec3& v, float angle, const vec3& axis)
{
	float new_x;
	float new_y;
	float new_z;

	float c, s, xx, xy, xz, yy, yz, zz, sx, sy, sz, ic;
	s = sinf(angle * DEG2RAD);
	c = cosf(angle * DEG2RAD);

	xx = axis.x * axis.x;
	xy = axis.x * axis.y;
	xz = axis.x * axis.z;
	yy = axis.y * axis.y;
	yz = axis.y * axis.z;
	zz = axis.z * axis.z;
	sx = s * axis.x;
	sy = s * axis.y;
	sz = s * axis.z;
	ic = 1.0f - c;

	new_x  = (xx * ic + c ) * v.x;
	new_x += (xy * ic - sz) * v.y;
	new_x += (xz * ic + sy) * v.z;

	new_y  = (xy * ic + sz) * v.x;
	new_y += (yy * ic + c ) * v.y;
	new_y += (yz * ic - sx) * v.z;

	new_z  = (xz * ic - sy) * v.x;
	new_z += (yz * ic + sx) * v.y;
	new_z += (zz * ic + c ) * v.z;

	v.x = new_x;
	v.y = new_y;
	v.z = new_z;

	v.normalize();
}

};	// namespace rtMath

#endif	// __MathLib_h__
